#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''
@Project ：V2 
@File    ：gm_fjsp_bp.py
@IDE     ：PyCharm 
@Author  ：郭星
@Date    ：2025/9/22 20:44 
'''
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
这份代码做了什么（简要）

两级+充装+QC：共 4 个阶段；每阶段有并机；QC 为组批机（容量/批时长）。

序列相关换型/吹扫：按“配方族”跳变加 setup_time 和 purge_gas（后者进入目标函数权重可扩展）。

安全相容（O₂/H₂）：若相邻禁忌，自动施加一段“等效隔离时间”（近似为插入N₂清洗段）。

有限缓冲：解码后估算缓冲轨迹，发现超额则把上游完工向后推（软阻塞迭代 2 轮，参数可调）。

有效期：以“充装完成→QC出库”的滞留时间与 shelf_life 比较；超出累加罚分。

能耗/电价：按阶段功率×时间×时段电价累加（中点取价）。

启发式/LNS：序列的交换/移位邻域，带“精英保留”（当前最好解）。

怎么扩展到论文版

严格缓冲/阻塞：把软阻塞改成离散时间桶（每 10min），加“库存平衡”不等式；或用 CP-SAT 的Reservoir/Cumulative。

一次→二次的物料平衡：现在把订单当子批直过全部工序；可拆为“先一次批（母气）→二次稀释批→订单充装子批”，并加物料守恒。

QC 更真实：多仪器、抽检比例、返工回流（在解码时生成返工任务再入队）。

目标多目标：切成 NSGA-II/多权情景；或者滚动优化加“变更惩罚”。

安全并发容量：用时间桶限制 
∑
𝑟
𝑟
𝑒
𝑐
𝑖
𝑝
𝑒
≤
𝑅
𝑚
𝑎
𝑥
∑r
recipe
	​

≤R
max
（跨工序合计），更贴近现场 EHS 管控。

如果你更偏向求最优（而非启发式），我可以把这套结构换成 CP-SAT（interval + sequence + reservoir + batch）的模型骨架；如果更偏向工业规模，我可以把 LNS 邻域做大（跨机迁移、簇交换、配方族整块重排），并加入多启动并行。
GM-FJSP-BP: Gas-Mix Flexible Job Shop with Buffers & Perishability
------------------------------------------------------------------
行业要素：
- 两级：一次混配(Mix1) -> 二次稀释(Dilute2) -> 充装(Fill) -> QC(组批)
- 序列相关换型/吹扫（按配方/浓度差）；禁忌相邻（如 O2/H2 必须N2隔离）
- 每工序后有限缓冲（满则前机阻塞：采用软阻塞迭代修正近似）
- QC 组批（容量上限、批处理时间）
- 有效期 & 分时电价
- 目标：加权（Cmax、迟期、吹扫耗气、能耗、电费、缓冲占用）

求解：
- 初解：按配方分组减少换型
- 解码：离散事件模拟 + 软阻塞迭代（修正前工序完工以消解缓冲超额）
- LNS邻域：交换/移位/跨机重分配，保留精英
"""

import math, random, numpy as np
from dataclasses import dataclass
from typing import List, Dict, Tuple, Optional
import matplotlib.pyplot as plt

# -------------------- 数据结构 --------------------

@dataclass
class Recipe:
    id: int
    name: str
    risk: float          # 安全并发权重（可用于容量约束，示例未激活）
    shelf_life: float    # 有效期（从充装完成算起）
    family: str          # 配方族（用于相容/分组）

@dataclass
class Order:
    id: int
    recipe_id: int
    qty: int             # 瓶数（本原型按“子批=订单”模拟；可扩展分拆）
    due: float           # 交期（小时）
    target: str          # 目标浓度标签（展示用途）

@dataclass
class StageConf:
    name: str
    m: int                  # 并机数
    rate: float             # （瓶/小时）或（单位量/小时）——此处用 qty/rate 得工时
    energy_kwh_per_h: float # 负荷功率
    buffer_cap: int         # 本工序后缓冲容量（件/子批）
    is_batch: bool=False    # QC为True：组批机
    batch_cap: int=0        # QC批容量上限
    batch_base: float=0.0   # QC基准时间
    batch_slope: float=0.0  # QC随样本数线性项（示例）

class Plant:
    def __init__(self, stages: List[StageConf], recipes: Dict[int, Recipe],
                 setup_time: Dict[Tuple[str, str], float],   # (stage_name, trans_key)->time
                 purge_gas: Dict[Tuple[str, str], float],    # 同上 -> 耗气
                 forbidden: Dict[str, List[Tuple[str,str]]], # stage_name -> [(from_family,to_family)]
                 price: Dict[int, float]):                   # 分时电价: hour->price
        self.stages = stages
        self.recipes = recipes
        self.setup_time = setup_time
        self.purge_gas = purge_gas
        self.forbidden = forbidden
        self.price = price

# 任务（子批）：以订单为粒度（可扩展为更细批）
@dataclass
class Lot:
    id: int
    order_id: int
    recipe_id: int
    qty: int
    due: float
    target: str

# 调度结果的片段
@dataclass
class Op:
    lot_id: int
    stage: int
    machine: int
    start: float
    end: float
    recipe_id: int
    label: str

# -------------------- 工具 --------------------

def seed_all(s=0):
    random.seed(s); np.random.seed(s)

def hour_price(price_map: Dict[int, float], t: float) -> float:
    h = int(math.floor(t)) % 24
    return price_map.get(h, price_map.get(-1, 1.0))

def plot_gantt(ops: List[Op], plant: Plant, title="Gantt"):
    if not ops: return
    by_stage = {}
    for op in ops:
        by_stage.setdefault((op.stage, op.machine), []).append(op)
    fig, axes = plt.subplots(len(plant.stages), 1, figsize=(12, 2.5*len(plant.stages)), sharex=True)
    if len(plant.stages)==1: axes=[axes]
    for s, stg in enumerate(plant.stages):
        ax = axes[s]
        for m in range(stg.m):
            key=(s,m)
            y = m
            ops_sm = sorted(by_stage.get(key, []), key=lambda x:x.start)
            for o in ops_sm:
                ax.barh(y, o.end-o.start, left=o.start, height=0.8)
                ax.text((o.start+o.end)/2, y, f"L{o.lot_id}", ha='center', va='center', fontsize=8, color='white')
        ax.set_yticks(range(stg.m)); ax.set_yticklabels([f"M{m}" for m in range(stg.m)])
        ax.set_title(f"{stg.name}")
        ax.grid(axis='x', alpha=0.3)
    axes[-1].set_xlabel("Time (h)")
    fig.suptitle(title)
    plt.tight_layout(); plt.show()

# -------------------- 解码（模拟） --------------------

class Decoder:
    """
    给定每个阶段的“序列”（lot_id 的顺序）和跨机分配策略（此处采用“随序就近分派”），
    进行多阶段离散事件模拟：
      - setup/吹扫：按相邻lot的“配方族”转移，叠加换型时间和耗气
      - 有限缓冲：初次解码后计算各阶段缓冲占用时序，如超过cap，软阻塞迭代：延后前一阶段完工
      - QC组批：按到达顺序组批，批处理时间 = base + slope * 样本数
      - 能耗电费：∫(功率 * 电价)·dt 近似为 Σ((end-start)*功率 * 电价(中点小时))
      - 有效期：充装完成到“假定出货时刻”的间隔超过 shelf_life 则算报废/处罚（演示为迟期罚）
    """
    def __init__(self, plant: Plant):
        self.plant = plant

    def setup_time(self, stage: int, from_rcp: Optional[int], to_rcp: int) -> Tuple[float, float]:
        if from_rcp is None: return 0.0, 0.0
        stg_name = self.plant.stages[stage].name
        f = self.plant.recipes[from_rcp].family
        t = self.plant.recipes[to_rcp].family
        key = (stg_name, f+"->"+t)
        return self.plant.setup_time.get(key, 0.0), self.plant.purge_gas.get(key, 0.0)

    def forbid_check(self, stage: int, from_rcp: Optional[int], to_rcp: int) -> bool:
        if from_rcp is None: return True
        stg_name = self.plant.stages[stage].name
        f = self.plant.recipes[from_rcp].family
        t = self.plant.recipes[to_rcp].family
        forb = self.plant.forbidden.get(stg_name, [])
        return (f,t) not in forb

    def proc_time(self, stage: int, lot: Lot) -> float:
        stg = self.plant.stages[stage]
        if stg.is_batch:
            return 0.0  # QC在组批阶段统一算时长
        # 简化：每子批时间 = qty / rate
        return lot.qty / stg.rate

    def decode_once(self, lots: Dict[int, Lot], seq_per_stage: List[List[int]]) -> Tuple[List[Op], Dict[int, float]]:
        """不考虑缓冲迭代，先按阶段独立排程（最早可开机），记录ops；返回ops和每lot在各阶段的完工时刻（仅用于QC等）"""
        ops: List[Op] = []
        lot_end_stage: Dict[Tuple[int,int], float] = {} # (lot_id, stage)->end
        last_recipe_on_machine = {}                    # (stage, m)->recipe_id

        # 逐阶段排
        for s, stg in enumerate(self.plant.stages):
            if stg.is_batch:
                # QC：按到达 Dilute2/Fill 完成后（上一阶段完工）顺序组批
                prev_s = s-1
                arrivals = [(lot_id, lot_end_stage[(lot_id, prev_s)]) for lot_id in seq_per_stage[s]]
                arrivals.sort(key=lambda x: x[1])
                # 单台或多台批处理机：简单轮转
                mac_ready = [0.0]*stg.m
                idx = 0
                while idx < len(arrivals):
                    # 取一批
                    take = arrivals[idx: idx+stg.batch_cap]
                    batch_start = max(min(mac_ready), max(a[1] for a in take))
                    batch_size = len(take)
                    batch_pt = stg.batch_base + stg.batch_slope*batch_size
                    m = int(np.argmin(mac_ready))
                    st = batch_start; en = st + batch_pt
                    for lot_id, _ in take:
                        ops.append(Op(lot_id, s, m, st, en, lots[lot_id].recipe_id, f"QC-{lots[lot_id].target}"))
                        lot_end_stage[(lot_id, s)] = en
                    mac_ready[m] = en
                    idx += stg.batch_cap
                continue

            # 非QC阶段
            mac_ready = [0.0]*stg.m
            prev_rcp_on_m = [None]*stg.m
            for lot_id in seq_per_stage[s]:
                lot = lots[lot_id]
                # 安全禁忌若不满足，则在序列中视作自动插入N2缓冲：等效增加额外setup（这里简化成大setup）
                # 选择最早能开机的机器（含setup），并满足安全
                best = None
                for m in range(stg.m):
                    # 预计基础开始：机器就绪 & 前序完工
                    prev_end = lot_end_stage.get((lot_id, s-1), 0.0)
                    st0 = max(mac_ready[m], prev_end)
                    # setup/吹扫
                    ok = self.forbid_check(s, prev_rcp_on_m[m], lot.recipe_id)
                    st_time, purge = self.setup_time(s, prev_rcp_on_m[m], lot.recipe_id)
                    # 若不安全，赋予很大的 setup 以“强制插入N2隔离”（或者你可以直接跳过该机）
                    if not ok: st_time += 1.0  # 1h 隔离等效
                    st = st0 + st_time
                    pt = self.proc_time(s, lot)
                    en = st + pt
                    if (best is None) or (en < best[0]):
                        best = (en, m, st0, st_time, purge)
                en, m, st0, st_time, purge = best
                st = st0 + st_time
                ops.append(Op(lot.id, s, m, st, en, lot.recipe_id, self.plant.stages[s].name))
                lot_end_stage[(lot.id, s)] = en
                mac_ready[m] = en
                prev_rcp_on_m[m] = lot.recipe_id

        # 记录每lot最后完工时刻
        lot_last_end = {}
        last_s = len(self.plant.stages)-1
        for lot_id in lots:
            lot_last_end[lot_id] = lot_end_stage[(lot_id, last_s)]
        return ops, lot_last_end

    def soft_blocking_iter(self, lots: Dict[int, Lot], seq_per_stage: List[List[int]], iters=2) -> List[Op]:
        """
        软阻塞迭代：先解码拿到ops，估算各缓冲时序；若缓冲超容量，
        就把“导致超额的上游工序完工”往后推（近似模拟阻塞），再重解码下游阶段。
        简化：最多做 iters 轮，小中规模下效果可接受；需要更严谨可用CP/离散时间桶。
        """
        ops, _ = self.decode_once(lots, seq_per_stage)

        for _ in range(iters):
            # 统计每阶段缓冲时间线：统计在[t, t+dt) 中的“已完成s但未开始s+1”的lot数
            # 近似：用每lot在 s 的 end_s 和 s+1 的 start_{s+1} （从ops中提取）
            by_lot_stage = {}
            for op in ops:
                by_lot_stage.setdefault((op.lot_id, op.stage), []).append(op)
            # 提取 start times 下游：构建映射
            start_next = {}
            end_curr = {}
            for op in ops:
                if op.stage < len(self.plant.stages)-1 and not self.plant.stages[op.stage+1].is_batch:
                    # 找下游同lot最早开始
                    nxt_ops = [o for o in ops if (o.lot_id==op.lot_id and o.stage==op.stage+1)]
                    if nxt_ops:
                        start_next[(op.lot_id, op.stage)] = min(o.start for o in nxt_ops)
                end_curr[(op.lot_id, op.stage)] = op.end

            # 对每阶段 s 的缓冲检查
            need_push = False
            push_tag = {}  # (lot_id, s)->delay
            for s, stg in enumerate(self.plant.stages[:-1]):  # 最后一阶段后无缓冲
                cap = stg.buffer_cap
                if cap <= 0: continue
                # 收集所有 (end_s, start_{s+1}) 区间
                intervals = []
                for lot_id in set([op.lot_id for op in ops if op.stage==s]):
                    if (lot_id, s) in end_curr and (lot_id, s) in start_next:
                        intervals.append((end_curr[(lot_id,s)], start_next[(lot_id,s)], lot_id))
                # 线扫（简化 dt）：统计任意时刻的同时滞留数是否>cap
                events=[]
                for a,b,lot_id in intervals:
                    if b<a: continue
                    events.append((a, +1, lot_id))
                    events.append((b, -1, lot_id))
                events.sort()
                cur=0
                for t,delta,lot_id in events:
                    cur += delta
                    if cur > cap:
                        # 缓冲溢出：把本lot在阶段 s 的完工推迟 Δ，使它的滞留晚点出现
                        push_tag[(lot_id, s)] = push_tag.get((lot_id, s), 0.0) + 0.2  # 推0.2h（参数可调）
                        need_push = True
                # end stage s
            if not need_push:
                break
            # 根据 push_tag 调整该stage的相应 lot 的“有效处理时长”（等价延迟完工）
            # 简化做法：在 decode 时为该 lot 在该 stage 额外加一个“延迟setup”
            # => 我们直接在 seq_per_stage 上不变，改成“临时setup矩阵”来重算
            # 这里为了不再创建复杂状态，采取二次解码并在完成时间后人为延后op.end
            new_ops=[]
            delayed = {(k[0],k[1]):d for k,d in push_tag.items()}
            for op in ops:
                if (op.lot_id, op.stage) in delayed:
                    d = delayed[(op.lot_id, op.stage)]
                    new_ops.append(Op(op.lot_id, op.stage, op.machine, op.start, op.end + d, op.recipe_id, op.label))
                else:
                    new_ops.append(op)
            ops = new_ops
            # 注：严格应级联影响下游start，这里下一轮 decode_once 会重新计算；为保持一致，我们重跑一次 decode_once
            ops, _ = self.decode_once(lots, seq_per_stage)

        return ops

    def kpis(self, ops: List[Op], lots: Dict[int, Lot]) -> Dict[str, float]:
        stg = self.plant.stages
        mks = max(o.end for o in ops) if ops else 0.0
        # 电费
        energy_cost=0.0
        for op in ops:
            pwr = stg[op.stage].energy_kwh_per_h
            mid = (op.start+op.end)/2
            energy_cost += (op.end-op.start)*pwr*hour_price(self.plant.price, mid)
        # 迟期
        tard=0.0
        by_lot = {}
        for op in ops:
            by_lot.setdefault(op.lot_id, 0.0)
            by_lot[op.lot_id] = max(by_lot[op.lot_id], op.end)
        for lot_id, fin in by_lot.items():
            due = lots[lot_id].due
            tard += max(0.0, fin-due)
        # 有效期（以充装完成->QC完成间隔代表滞留；过期计罚）
        per_pen=0.0
        # 取每lot在充装(stage=2)完工 与 最终出库(=QC完工 stage=3)的差
        for lot_id in lots:
            fill_end = max([o.end for o in ops if o.lot_id==lot_id and o.stage==2], default=0.0)
            ship_end = max([o.end for o in ops if o.lot_id==lot_id and o.stage==3], default=fill_end)
            life = self.plant.recipes[lots[lot_id].recipe_id].shelf_life
            over = max(0.0, ship_end - fill_end - life)
            per_pen += over
        return dict(Cmax=mks, EnergyCost=energy_cost, Tardiness=tard, PerishPenalty=per_pen)

# -------------------- 解编码（解的表示与邻域） --------------------

class Solution:
    def __init__(self, seq_per_stage: List[List[int]]):
        self.seq = [list(s) for s in seq_per_stage]
        self.obj = None
        self.ops: List[Op] = []

class Heuristic:
    def __init__(self, plant: Plant, lots: Dict[int, Lot], decoder: Decoder):
        self.plant=plant; self.lots=lots; self.dec=decoder

    def initial(self) -> Solution:
        # 初解：各阶段按配方族分组（减少换型），同族内部按交期
        by_family = {}
        for lot in self.lots.values():
            fam = self.plant.recipes[lot.recipe_id].family
            by_family.setdefault(fam, []).append(lot)
        for fam in by_family:
            by_family[fam].sort(key=lambda x:x.due)
        init_seq=[]
        for s,_ in enumerate(self.plant.stages):
            seq=[]
            for fam in sorted(by_family.keys()):
                seq += [lot.id for lot in by_family[fam]]
            init_seq.append(seq)
        sol = Solution(init_seq)
        sol = self.evaluate(sol)
        return sol

    def evaluate(self, sol: Solution) -> Solution:
        ops = self.dec.soft_blocking_iter(self.lots, sol.seq, iters=2)
        k = self.dec.kpis(ops, self.lots)
        # 目标加权
        w = dict(Cmax=1.0, Tardiness=2.0, EnergyCost=0.01, PerishPenalty=5.0)
        obj = sum(w[kx]*k[kx] for kx in w)
        sol.obj = obj
        sol.ops = ops
        return sol

    # 邻域：在随机阶段做交换/移位；或在两台机之间做交叉（本原型序列级，不显式机器层）
    def neighbor(self, sol: Solution) -> Solution:
        new = Solution(sol.seq)
        s = random.randrange(len(new.seq))
        if len(new.seq[s])>=2:
            if random.random()<0.5:
                i,j = sorted(random.sample(range(len(new.seq[s])),2))
                new.seq[s][i], new.seq[s][j] = new.seq[s][j], new.seq[s][i]
            else:
                i = random.randrange(len(new.seq[s]))
                j = random.randrange(len(new.seq[s]))
                x = new.seq[s].pop(i); new.seq[s].insert(j, x)
        return new

    def lns(self, iters=5000, temp0=0.1) -> Solution:
        cur = self.initial()
        best = cur
        for t in range(1, iters+1):
            cand = self.neighbor(cur)
            cand = self.evaluate(cand)
            if cand.obj <= cur.obj or random.random()<math.exp((cur.obj-cand.obj)/max(1e-9, temp0)):
                cur = cand
                if cand.obj < best.obj:
                    best = cand
            if t%500==0:
                print(f"[LNS] iter {t}  best={best.obj:.2f}")
        return best

# -------------------- 数据生成（可复现实验） --------------------

def build_demo(seed=1) -> Tuple[Plant, Dict[int, Lot]]:
    seed_all(seed)
    # 配方（示例：O2、H2、N2 family 及若干目标）
    recipes = {
        0: Recipe(0,"O2-40%", risk=2.0, shelf_life=30.0, family="O2"),
        1: Recipe(1,"O2-2%",  risk=1.5, shelf_life=20.0, family="O2"),
        2: Recipe(2,"H2-5%",  risk=3.0, shelf_life=15.0, family="H2"),
        3: Recipe(3,"N2",     risk=0.5, shelf_life=60.0, family="N2"),
        4: Recipe(4,"CO2-10%",risk=1.0, shelf_life=40.0, family="CO2"),
    }
    # 阶段：Mix1, Dilute2, Fill, QC
    stages = [
        StageConf("Mix1",    m=2, rate=200.0, energy_kwh_per_h=5.0, buffer_cap=4),
        StageConf("Dilute2", m=2, rate=150.0, energy_kwh_per_h=4.0, buffer_cap=5),
        StageConf("Fill",    m=3, rate=100.0, energy_kwh_per_h=3.0, buffer_cap=8),
        StageConf("QC",      m=1, rate=0.0,   energy_kwh_per_h=2.0, buffer_cap=999,
                  is_batch=True, batch_cap=5, batch_base=0.3, batch_slope=0.08),
    ]
    # setup/吹扫（按配方族跳变），仅示例性设置
    setup_time={}
    purge_gas={}
    families = ["O2","H2","N2","CO2"]
    for s in ["Mix1","Dilute2","Fill"]:
        for f in families:
            for t in families:
                delta = 0.05 if f==t else (0.2 if {f,t}=={"O2","N2"} else 0.4)
                setup_time[(s, f+"->"+t)] = delta  # 小时
                purge_gas[(s, f+"->"+t)] = 0.2 if f!=t else 0.05
    # 禁忌（需要N2隔离）：O2<->H2 在所有工序
    forbidden = { s.name: [("O2","H2"), ("H2","O2")] for s in stages }
    # 分时电价（简单峰平谷）
    price = {h:(1.2 if 9<=h<12 or 18<=h<22 else (0.8 if 0<=h<6 else 1.0)) for h in range(24)}

    plant = Plant(stages, recipes, setup_time, purge_gas, forbidden, price)

    # 订单 -> lot（每订单一个子批；可扩展为订单拆分）
    lots={}
    O=18
    for i in range(O):
        r = np.random.choice(list(recipes.keys()))
        qty = np.random.randint(80, 220)    # 瓶
        due = np.random.uniform(10, 50)     # 小时
        target = recipes[r].name
        lots[i] = Lot(i, i, r, qty, due, target)
    return plant, lots

# -------------------- 主程序 --------------------

def main():
    plant, lots = build_demo(seed=7)
    dec = Decoder(plant)
    heu = Heuristic(plant, lots, dec)
    best = heu.lns(iters=3000, temp0=0.05)
    print("\n== KPIs ==")
    k = dec.kpis(best.ops, lots)
    for kx,v in k.items():
        print(f"{kx:>14s}: {v:.3f}")
    print(f"{'Objective':>14s}: {best.obj:.3f}")
    plot_gantt(best.ops, plant, title="GM-FJSP-BP (Heuristic)")

if __name__ == "__main__":
    main()
